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We perform a numerical simulation of Faraday waves forced with two-frequency oscillations using 
a level-set method with Lagrangian-particle corrections (particle level-set method). After validat¬ 
ing the simulation with the linear stability analysis, we show that square, hexagonal and rhom- 
boidal patterns are reproduced in agreement with the laboratory experiments [Arbell and Fineberg, 
Phys. Rev. Lett. 84 , 654 (2000) and Phys. Rev. Lett. 85 , 756 (2000)]. We also show that the 
particle level-set’s high degree of conservation of volume is necessary in the simulations. The nu¬ 
merical results of the rhomboidal states are compared with weakly nonlinear analysis. Difficulty in 
simulating other patterns of the two-frequency forced Faraday waves is discussed. 


I. INTRODUCTION 

Faraday waves^, known to exhibit various kinds of crystalline patterns in simple settings, have attracted 
many researchers for about two hundred years. Faraday waves are the surface waves between two superposed 
immiscible fluid layers subjected to a vertical vibration. Even recently astounding exotic phenomena continue 
to be found in laboratory experiments on Faraday waves. For example, in Faraday waves with a certain non- 
Newtonian fluid (shear-thickening fluid), the behavior of the interface is far beyond what one can imagine 
from the interface motion between air and water—. In another surprising experiment, a droplet slightly 
submerged in a liquid substrate under a vertical oscillation is found to behave dynamically like a snaked. 
To physically understand these phenomena, numerical simulations of them, which may not be possible now, 
are expected to play a decisive role. 

As a first step to build such numerical methods, we here study numerically Faraday waves subjected to 
a two-frequency forcing in a Newtonian fluid. There are a number of experimental results with this forcing 
setting^^— , where much richer variations of the selected patterns are found than in the single-frequency 
forced cases as listed below. 

The study of two-frequency forced Faraday waves starts with the experiments by Edwards et and 
Muller—. The two-frequency forcing can be written as Ai cos(mujot) A 2 cos(ncoot -h (f>) and characterized by 
the integers m and n. Edwards et al. explored various ratios of the two frequencies such as m : n = 3 : 5,4 : 
5, 6 : 7,4 : 7 and 8 : 9 and mainly investigated the ratio 4 : 5. They observed the quasi pattern, which has 
a long-range orientational order but no spatial periodicity. On the other hand, the experiment by Muller is 
focused on the driving ratio of 1 : 2 and produces a triangular pattern. 

In the linear regime of the two-frequency forced case a bicritical point exists at which two normal modes 
with different wavenumber moduli become simultaneously unstable (for the single-frequency forced case the 
bicritical point can be formed by tuning the frequency of the forcing for shallow layers^). The unstable 
modes then interact with each other nonlinearly. In the neighborhood of the bicritical point many complex 
patterns are expected to be found. A number of experiments around the bicritical point were conducted by 
Kudrolli et al.—, Arbell et al^A^L and Epstein et al.— Kudrolli et al. observed patterns that they named 
superlattice-1 and superlattice-2. Arbell et al. and Epstein et al. observed double hexagonal superlattice 
(DHS), subharmonic superlattice states (SSS), oscillon, two-mode superlattices (2MS) and 2k rhomboidal 
states (2fcR). Each pattern can be characterized by the number of excited (discrete) Eourier modes and by 
the nonlinear resonance among them. 

To the best of our knowledge, numerical simulation of the two-frequency forced Faraday waves based 
on the Navier-Stokes equations solving the motion of both the top and bottom fluids is reported in this 
paper for the first time. However, such a simulation, not limited to the two-frequency forced case, requires 
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treatment of the interface with the surface tension force. We therefore must employ one of the interface¬ 
tracking schemes such as the volume-of-fluid methods, the level-set methods, the front-tracking methods, 
(see, e.g., an advanced textboobi^). In this study, we adopt the level-set method and investigate whether or 
not the simulation of two-frequency forced Faraday waves with the level-set method is consistent with the 
experimental results. The reason for adopting the level-set method will be described later. 

The first numerical simulation of the single-frequency forced Faraday waves in three dimensions was 
performed by Perinet et al.—, who reproduced the square and hexagonal patterns in quantitative agreement 
with the laboratory experiment by Kityk et aL— . Perinet et al.— used a front-tracking method. It is 
necessary, for example in simulating oscillon or snake-like patterns, to allow for overturning and topological 
change of the interface. We hence believe that other interface-tracking schemes should be explored and 
tested for a wider class of the Faraday waves. Another numerical issue concerns the density difference 
between the top and the bottom fluids. In typical laboratory experiments, these are air and water at room 
temperature, meaning three orders of magnitude difference in the densities. To handle this large difference, 
it is known that a high quality solver for the pressure Poisson equation is needed regardless of the choice of 
interface-tracking scheme^. We use a preconditioned BiCGSTAB. 

On the theoretical front of the two-frequency forced Faraday waves, linear stability analysis and weakly 
nonlinear theory are available. Linear analysis was performed by Besson et al.—, which is an extension of 
the single-frequency forced casei^,. Their resultsi^ agree with the experiments quantitatively. In the weakly 
nonlinear analysis, whose emphasis is on the pattern selection of the two-frequency Faraday waves, Silber 
et al, Tse et al., Porter et al. and Topaz et al—^— formulated an amplitude equation up to third order in 
amplitude by applying symmetry based arguments. 

By analyzing the structure of the three-wave resonance, they succeeded in explaining many selected 
patterns qualitatively. Quantitative prediction of the pattern can be obtained if the amplitude equation of 
the two-frequency forced Faraday waves is derived from the Navier-Stokes equation with a realistic boundary 
condition. However this is a formidable task. A reduced hydrodynamic equation of the two-frequency 
Faraday waves was derived by Zhang et al—. From this reduced equation, the amplitude equations are 
derived and analyzed by assuming infinite depth and small viscosity^^— . Weakly nonlinear analysis based 
on the Navier-Stokes equations with infinite depth was carried out by Skeldon et al—. This approach with 
realistic amplitude equations is successful in explaining many patterns observed in the two-frequency forced 
Faraday waves. Nevertheless, there are some patterns, such as oscillonsiS, which are not explained so far 
by the weakly nonlinear analysis. In the effort to understand these patterns, numerical simulation of the 
Faraday waves plays a complementary role. 

For this reason, we develop a method of numerical simulation of the two-frequency forced Faraday waves, 
which is consistent with the experiments. Specifically, we here simulate three patterns observed in the 
experiments by Arbell et a/.— In particular the rhomboidal pattern does not appear in the single-frequency 
forced Faraday waves. In order to validate the simulations, we compare our results with the linear stability 
analysis of two frequency Faraday waves^^. Next, in the nonlinear regime, we reproduce the square pattern 
and the hexagonal pattern with the same physical parameters as the respective experiments. After that, 
we reproduce and study the rhomboidal state. During the simulations, we compare two kinds of level- 
set methods: one is the original implementatio n^^i^° and the other is the level-set method with Lagrangian 
particles (particle level-set method)^. Finally, we discuss the difficulty of simulating other patterns observed 
in the experiments. 

The organization of the paper is the following. In Section [III we describe the fluid dynamical equations of 
the Faraday waves, the two level-set methods and numerical discretization of the equations. The numerical 
results are presented in Section fllll More specifically, comparisons of the simulation with the linear analysis 
and simple patterns such as square and hexagonal patterns are presented in Section IIII Al and IIIIB1 The 
simulation of the rhomboidal states is shown in Section nirn In Section IIIIDl we compare the original 
level-set method and the particle level-set method. Our summary and discussion are in Section llVI 


II. EQUATIONS AND NUMERICAL METHOD 


In this section, we describe our numerical method for the governing equations and the boundary conditions 
used for simulating Faraday waves oscillated by the two-frequency forcing. 
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A. Navier-Stokes equations 


Faraday waves occur on the interface between an upper and a lower immiscible fluids. We employ the 
one-fluid description of the problem. Numerically we simulate the dynamics in both fluid layers. The 
incompressible Navier-Stokes equations are written as 

V - M = 0, (1) 

pDtu = — Vp -I- pG -I- V • p(Vm + Vtt^) -I- s. (2) 

Here, Dt^p, u are the material derivative, the pressure and the velocity, and s, p and p are the surface force, 
the density and the viscosity, respectively. The vector G is the gravitational term in the reference frame of 
the container. 


G = {—g -\- Ai cos(mwot) -I- A 2 cos{nujot + 0))ez (3) 

where g, Ai, H 2 , wq, 0, are the gravitational acceleration, the amplitude of the first periodic forcing, 
the amplitude of the second periodic forcing, the base angular frequency of the periodic forcing, the phase 
shift between the two modes, the unit vector in the vertical z-direction. In this paper, we set the integers 
m, n to 2, 3. We also use the notations = muiQ = 2uJo, UJ 2 = nujQ = 3aJo- 

On the top and bottom boundaries, no-slip boundary conditions are assumed. For the horizontal direction, 
we assume periodic boundary conditions. The interface location z = C,{x, y, t) obeys the kinematic boundary 
condition. In term of this C,, the density p and p are written as: 




[pt.pt) z > C(a:,2/,t), 
{Pb,Pb) z<C{x,y,t), 


( 4 ) 


where pt, pt are the density and the viscosity of the top fluid and pb are the density and the viscosity 
of the bottom fluid. 

In this sharp interface description, the density and the viscosity change discontinuously at the dynamically 
evolving interface. This situation is a challenge for numerical simulations. To circumvent this difficulty, 
various numerical methods have been proposed, such as the volume-of-fluid methods, the level-set methods 
and the front-tracking methods, just to name a fe w^^i^^ . In this study, we adopt the level-set method. The 
reason is as follows. The level-set method has a high numerical accuracy of the normal vector and the 
curvature of the interface, hence adequate for the gravity-capillary waves. However, it is well known that 
the level-set method does not have good mass conservation properties^. A number of improvements have 
been proposed^^— . Among them, we use the level-set method corrected with Lagrangian particles, the 
so-called particle level-set method, to ensure volume conservation^. This conservation problem is discussed 
in detail in Section PlI PI 

In the following section ITl B[ we describe the level-set method without the particles, here we call the original 
level-set method, and the particle level-set method and their numerical discretizations. The description of 
the discretization of the Navier-Stokes equations follows later. 


B. level-set method 

1. level-set function 

We use the level-set approach^^ to describe the interface motion. Here the level-set function the 

signed distance from the interface, indicates the interface. We define </> > 0 in the top fluid and ^ < 0 in the 
bottom fluid. The level-set function obeys the following equation 

dt^+{u-V)(l) = 0, (5) 

which is discretized with the 5th-order WENO scheme^! and integrated in time through the 3rd-order TVD 
Rung-Kutta method^. 
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2. Reinitialization of level-set function 

It is known that the analytic integration of Eq. dSD does not ensure that is the signed distance 

function from the interface. By definition, being the signed distance function requires |V(/)| = 1. However 
this unit gradient condition is not satisfied since the Lagrange derivative of |V(/)p is not zero but = 

—2(V0)(Vm)(V0). To enforce the condition (in practice, we do so just around the interface), the distance 
function is re-initialized at each time step from the following initial value problem with the virtual time r— 

r)fj 

— =sgn((/))(l-|Vd|)-f A/(())), (6) 

d{x,y,z,T = 0) = (j){x,y,z). 

Although we call r virtual time, its dimension is length. Ideally, the function d{x,T) as r —>■ oo gives the 
corrected signed distance function for all the computational domain. Here, we set 4’{x,t) = d{x,T = ti) for 
some value t;. This t; corresponds to the largest distance from the interface to which we demand (j) be the 
signed distance. In this paper, we use ti = e, where e is the half width of the diffuse interface and set to 
2Az, where Az is the grid spacing in the vertical z-direction. The functions A(a:), /(</>) in Eq. ([5]) are given 
as 


A(a;) = - 

H'{(j))L{(j), d)dx 


(7) 

Ini.)H'i^)fi^)dx ’ 


/(</>) = id'(</>)! Vc^l, 


(8) 

" = If 


(9) 

1 

fo 

if (/) < — e. 


ii(<^) = 

i{l + | + isin(^)} 

if — e < 0 < e, 

(10) 

1 


if (/) > e, 


where H(a:) is a small region centered at the point x, L{(p,d) = sgn((())(l — |V(i|), e = 2Az is the prescribed 
interface width, H{(j)) is the smoothed Heaviside function and H'{(j)) is the smoothed delta function. 


Numerically, the reinitialization is done in the following way. Firstly, we ignore the term X{x)f((f>) in the 
Eq. ([6|) and solve 


— =sgn((/))(l-|Vd|), (11) 

where the discretization in space is the same as that of Eq. The integration in the virtual time is 
discretized as follows. 


d' = d" -f ArL(((>,d”) 

d* = i(3d”-bd' +ArL((/),d')) 

d^n+i = i(d" + 2{d* + AtL{ 4>, d*))), 

where At = 0.5 min(Aa;, Ay, Az). 

Secondly, we calculate A(a;) according to the following equation 


( 12 ) 

(13) 

(14) 


(15) 


Here Xijk denotes X{xijk) on the grid point Xijk specified by the index (i,j,k). The integral range flijk 
describes the cell region associated with the grid point. For the three dimensional case, by following the 
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two-dimensional version^, we discretize the integral of some function g(x) in the cell as 


f /\t./\ y/\ 7. 

/ g{x) dx = —-b 


Ij'+ljfc “1“ 1,^ — 1,A; “1“ 

+5z+l,j,fc+l + 

+5i-lj,fe-l) + 88((?i+ij_fc + gi-l,j,k 
+gi,j+l,k + gi,j-l,k + gi,j,k+l 

+gi,j,k-i) + 876gij^k], 

where Ax, Ay, Az are the grid spacings along the x, y, z directions. 

Finally, is calculated by 

d"+i = d'"+i -b AT\ijkH\(l))\yct)\. 

In practice, we take the total number of the virtual time steps as tj/At ~ 4. 


(16) 


(17) 


C. Particle level-set method 

In order to improve the volume conservation of the level-set method, it has been proposed to utilize 
Lagrangian information to correct the level-set function by adding marker particles near the interface. Our 
procedure of the particle level-set method is basically the same as that of Enright et al^. The differences 
are in the error correction and the reseeding strategy. 


1. Initialization of particles 

The marker particles are spread in the neighborhood of the interface, in which \(j)\ < 3max(Aa;, Ay, Az) 
is satisfied. The number of particles in each cell is set to 64. A marker particle has sign Sp = 1 or —1 
and the radius Vp. There are a number of strategies for setting the sign and radius. One simple strategy is 
to set the sign to that of the level-set function at the particle and the radius to the absolute value of the 
level-set function. However, we follow the more sophisticated strategy proposed by Enright et al. to improve 
numerical results. 

The strategy is as follows. Initially, the particle’s sign is set randomly. In order to have the same sign 
between the particle and the level-set function, the particles at Xp, (p = 1, 2, 3,...) are iteratively moved by 
the following recurrence relation 

x;+^ = x; + 2--{<kpoai - Hxn)N{x;), (is) 

where N{Xp) = is the normal vector. Here (j)goai is set as follows. The sign, sgn((^p), is set to have 

the same as that of the particle Sp. In addition, the absolute value is chosen to be a uniformly distributed 
random variable in the range &i„in < \(l>goai\ < &max- In this study, 6min is set to 0.1 min(Aa:, Ay, Ax) 
and 6max is set to 3 max(Ax, Ay, Az). Each particle is moved repeatedly by Eq. (fT^ until it satisfies the 
condition &i„in < Sp4>(xp) < ^max- Finally, each particle radius is set according to 

{ Tu SpCfyiXp) > Tu, 

Sp(l}{Xp) Tl Spij^iXp) < Tu, (il^) 

n Sp(t){Xp) < ri, 

where ri and are lower and upper limits of particle radius to prevent the creation of particles which are 
too small or too large. We use n = 0.1 min(Aa;, Ay, Az) and = 5n. The particle radius Xp is used to 
correct the level-set function later. After this procedure, the positive particles at position Xp are in the 
4>{xp) > 0 side (the top fluid) and the negative particles are in the (j){xp) < 0 side (the bottom fluid). The 
envelope formed by the circles of the same-sign particles coincides with the interface (f> = 0. 
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2. Advection of particles 

Each particle at position Xp{t) is advected by 

( 20 ) 

The velocity at the particle position u(xp(t),t) is calculated with trilinear interpolation from the velocity 
vectors on the nearby cell faces. The 3rd order TVD Runge-Kutta method is used to integrate Eq. dsni) in 
time. 


3. Error correction of level-set function 


As a result of the advection Eq. (uni, some particles move across the interface ^ = 0. Such escaped 
particles are used to correct the level-set function cj) in the following manner. First, particles placed on the 
wrong side {(t>{xp) x Sp < 0) are considered to have escaped. Second, we introduce the signed distance 
function between the escaped particle and a point x, which is calculated with the particle radius Vp as 

(j)p{x) = Sp{rp - \x - Xp\). (21) 

This signed distance is positive {(j)p{x) > 0) if the point x is within the positive ball {sp = 1) of radius Vp 
centered on Xp. The distance function corrected by the escaped positive (negative) particle ())+(£c) {4>-{x)) 
is calculated from 


(j)+{x) = max {(j)p,(j){x)), 
(j)-{x) = min {(j)p,(j){x)). 

Wp^E~ 


( 22 ) 

(23) 


Here and E denote the sets of the escaped positive and negative particles. Finally, the level-set function 
(j){x) is corrected as 


cj){x) 


(j)+{x) if |(;i+(a;)| < |(;i_(a;)|, 
(/)_(a;) if |(;i_(a;)| < |(/)+(a;)|. 


(24) 


Ideally, after this correction of the distance function, all the particles tagged as escaped have the same sign 
as the corrected distance function. However, with our implementation of the correction in the preliminary 
calculations, we find that some particles do not have the same sign of the corrected distance function. If we 
use such particles with the wrong sign in the next correction process, the interface becomes nearly singular, 
which we consider a numerical artifact. Therefore we ignore such escaped particles in the later correction 
processes. The point differs from the usual procedure of the particle level-set method^. 


4. Reseeding of particles 


Generally, as a result of advection of the particles by a flow, some regions lack sufficient particles to 
correct the level-set function. We reseed the particles where needed. Specifically, in the cells near the 
interface {\(j){x)\ < bmax), we keep the number of particles in a cell to 64 by adding particles for cells which 
particles exit or deleting particles for cells which particles enter. In our simulation, the reseeding procedure 
is executed with the following two strategies. The first strategy is that the reseeding is done after 40 time 
steps from the previous reseeding. The second strategy is that the reseeding is done when the surface area 
of the interface increases by 30% after the previous reseeding time. The surface area is calculated from 


A = j 6 {<j))\V(l)\dx, 


(25) 


5W = 



(l-bCOs(^)) 


kl > e, 

l<^l < e, 


(26) 


where the smoothed delta function is the same as of Eq. ([9]). 

In summary, the one-step update of the level-set function with particles is carried out by the following 
steps^. 
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1. The level-set function is advected by Eq. 

2. The particles are advected by Eq. dini). 

3. The error of the level-set function is corrected by the procedure described in the Section III C 31 

4. The level-set function is reinitialized as described in Section Pi B 21 

5. The error of the level-set function is once more corrected by the particles as described in Section Pi C 31 
For the original level-set method without particles, the second, third and fifth processes are omitted. 


D. Discretization of Navier-Stokes Equations 


We use the following temporal discretization of the incompressible Navier-Stokes Eqs. o and with 
the projection method and with adaptive time stepping 

u* =it" -I- At' 

+ ( 1 


1-b 

At®®-! ^ 


Ay-®®-! 

n-1 

At®® 2 

jA + 

At®® 

1 + 

At®®-^ ^ 

)d:- 

Ay®®-! 

nn-! 

At®® , 

At®® 


=u* 




pU 




(27) 

(28) 


Here the superscript n denotes the value at the n-th time step t" = which At® is the time step 

size for the t-th step. We discuss later how to determine them. A" is the advective term, Z)" and D* are 
the viscous terms involving the same component as on the left-hand-side and at the intermediate step, D" 
is the viscous term involving the other components, G" is gravity and is the surface force term. Here, 
by * we denote the intermediate step. As in the standard way of the projection method, the pressure term 
is calculated from the divergence free condition. Equation (l28l) acted upon by V- becomes 

(29) 


- ^ = V f—Vp"+i 
At®* Vp” 

The advection term A" is described by 

A" = m" • Vm". 

The x-component of the viscous term Z>" to be treated implicitly is 


= 


1 


P 


V (ry"V<) + ^ I ry 


A 

dx 


dul 

dx 


(30) 


(31) 


The viscous term D*^ is defined similarly but with the intermediate velocity u*. The x-component of the 
viscous term G" to be treated explicitly is 


£)" = — 
ex „ 


d 


.du^ 

dx 


d 


dx 


(32) 


Other components of the viscous terms are described in the same manner. The gravitational term G" is 

G" = {—g + Ai cos(a;it") -I- A 2 cos(a; 2 t" + 0))e.z. (33) 

The surface force term S’® is 


S" = 




K®® = V - n®®. 


(34) 
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Regarding the spacial discretizations, the advection term, A", is discretized with the 2nd-order ENO scheme. 
The other derivative terms in £)*,£>", ZD", n”, k" are discretized with the 2nd-order central difference 
scheme. 

Concerning the boundary conditions, we assume periodic boundary conditions in the horizontal directions 
{x and y directions). For the vertical z direction, we assume the non-slip condition at z = 0,2^2 


w|z=o,L^ = 0- 

The boundary condition for the pressure in solving the Poisson equation (1291) is 


(35) 


dp 

dz 


= 0 . 

2: —0, Lz 


(36) 


For calculating u* in Eq. (|27|. the localized ILU preconditioned BiCGSTAB method^ is adopted and the 
Poisson equation of the pressure is solved by the multigrid preconditioned BiCGSTAB method^^. 

Finally, we describe how we determine the variable time step size At* which is determined by 


Af = cmin(At 5 , Aty{x), Atrj{x), Atlfi{x)), 


(37) 


where we set the safety constant c = 0.40. Flere Ats, At/ and At,, are the time scales of surface force, the 
vertical vibration and the viscosity, respectively. The time scale Atcfi concerns the CFL condition. These 
reference time scales are defined as 


Ats = 


{pt + Pb)dh^ 


4707 


, (dti = min(Aa;, Ay, Az)), 


Affix) -^Qi , ^ 1 


1 


1 + 1 ’ 
Ax-‘ ' ' Az‘‘ 


"x I **1/ I K 
Ax ^ Ay Az 


Typically Ats is the smallest in our all simulations. 


(38) 

(39) 

(40) 


III. NUMERICAL RESULTS 

Our goal in this paper is numerical simulation of the 2k rhomboidal states observed in the laboratory 
experiments by Arbell et alM^. To the best of our knowledge, the rhomboidal states have not previously 
been obtained in numerical simulations of the Navier-Stokes equations. In particular, we use the same 
bulk fluid parameters as the experiments. The only difference is the geometry of the system, i.e., the 
domain size and the boundary conditions. In the simulations we apply periodic boundary conditions in the 
horizontal directions, with which we can reduce numerical cost by not simulating many repeated patterns 
in the computational domain. While the experiments are conducted in an open container, we assume the 
presence of a rigid wall above the top fluid, on which the no-slip boundary condition is applied. This setting 
is numerically easier than simulating the top-fluid motion in a semi-infinite domain. We assume that the 
top-fluid height is four times larger than the bottom-fluid depth. 

Before presenting the simulation, we first describe the validation of our fully nonlinear simulation with the 
above mentioned geometry by comparing with the linear stability analysis of two-frequency forced Faraday 
wavesi^. This validation process is the same as Perinet et The second test then is to reproduce 

the square and hexagonal patterns observed in the same two-frequency forced experiments^di. A direct 
numerical simulation of the square and hexagon patterns for the single-frequency forced Faraday waves is 
performed by Perinet et The result on the rhomboidal states is presented after the validations. 


A. Comparison with linear stability analysis 

We now compare the critical amplitudes of the oscillations calculated with the fully nonlinear numerical 
simulation with those calculated with the linear stability analysis. We write the two-frequency forcing 
as a(cos(x) cos(wit) -|- sin(y) cos(w 2 ^ + d)). In the linear analysis, once we fix the physical parameters as 
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Table I. Parameter values for the linear stability analysis. These parameters (except for and 9) are identical to 
the experiment by Arbell et al^ 


pt 

1.293 

[kgm-®] 

Pb 

9.500 X 10 ^ 

[kgm-3] 

Vt 

1.822 X 10 "® 

[kgm“^ s“^] 

Vb 

2.185 X 10 "^ 

[kgm-W-i] 

oji = 2 luo 

3.770 X 10 ^ 

[s-^] 

UJ2 = ScjQ 

5.655 X 10 ^ 

[s-i] 

6 

0 

[rad] 

a 

2.150 X 10 "^ 

[kgm-i] 

g 

9.807 

[ms-2] 

L, 

1.00 X 10 "^ 

[m] 

bottom-fluid depth 

2.00 X 10 "® 

[m] 


shown in Table |I] and the mixing angle y, then the critical value of a, denoted as Uc, and the associated 
critical wave number can be calculate d^^d^ . Here we assume that either harmonic frequency (uji,uj 2 ) or 
sub-harmonic frequency (wi/2, W 2 / 2 ) gives the lowest critical value. The critical amplitudes Oic = acCOs(x) 
and a 2 c = Oc sin(x) for the mixing angle y from 0° to 90° are shown as the solid line in Fig. [TJ 

Meanwhile, with the particle level-set simulation, we determine the critical amplitudes by adding small 
perturbations to basic modes for ten different values of the mixing angle, y = 0°, 10°, 20°,..., 90°. The 
results are denoted as points in Fig. [TJ The way to estimate the critical amplitudes aic and a 2 c in the 
nonlinear simulation is as follows: (i) the perturbation is added to the normal mode whose wavenumber is 
set to the critical wave number (fee) calculated from the linear analysis. More precisely, the perturbation, 
^sin(fcca::), is added to the flat interface, where the perturbation amplitude ^ is set to 2 x 10“^Lz. (ii) the 
interface height ^(x,y,t) at the center of the calculation domain is monitored throughout the simulation 
for given a and y. The interface height z = is calculated from the zero points of the level-set 

function (j){x,y, z^t) = 0. We perform such simulations by changing a and estimate the critical value Qc- 
More precisely, we take the absolute relative difference between the two peak interface heights at t = 4.44r„ 
and 6.44T„ where T„ = 2ttIujq = 27r/(a;i/2) = 27r/(a;2/3). Then, 2r„ is the minimal period of the two waves 
with the subharmonic frequencies uJi/2 and W 2 / 2 . If the difference is smaller than 10“^, this a is regarded 
as the critical amplitude of the level-set numerical simulation. 

As shown in Fig. [1] the critical amplitudes calculated with the level-set simulation tend to be greater 
than those calculated with the linear analysis in the 02 dominant region, namely ai/g < 2.6. The absolute 
relative error between the linear analysis (line) and the simulation (point) in the region is about 0.02. The 
agreement between the two results is hence obtained with two-digit accuracy. 


B. Square and hexagonal patterns 

Having validated the simulation of the two-frequency forced Faraday waves in the linear regime, we now 
move to two nonlinear cases: the square and hexagonal patterns. Note that the two patterns are also 
observed in the single-frequency Faraday waves. 

First we reproduce the square pattern observed in the experiments by Arbell et The physical 

parameters are shown in Table El We select the amplitudes of the forcing Ai and A 2 according to the 
following reasons: (i) we aim to conduct the simulation in the weakly nonlinear regime; (ii) we aim to 
set the values to be neither close to nor far from the bicritical point. The selected values of Ai and A 2 
in Table El are of course in the square-pattern domain of the phase diagram obtained experimentallyi^. 
However the geometry of the simulation is different. We set the lateral dimensions of the computational 
domain so that it includes one square 27r/fci = Ai, where fci is the critical wave number, which is found to be 
1.436 X 10^ from the linear stability analysis described in the previous section. In other words, we set 

the computational domain to a square box with = Ly = \i. This setting is the minimal computational 
domain which supports the periodic square pattern. The number of grid points used in each horizontal 
direction is denoted by = Ny. The square pattern consists of the four discrete Fourier modes shown as 
black dots in Fig.jSKa). These modes, called resonant modes, are on the circle of radius ki. The amplitude 
of the resonant wavevectors can be calculated from the linear stability analysis. The direction of those can 
be estimated from the experimental data. The experimental information of the direction is trivial in the 
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Figure 1. Comparison of the critical vibration amplitudes obtained with the numerical simulation to those calculated 
with the linear stability analysis. 


Table II. Parameter values for the square pattern. The number of the grid points for each direction are N^, Ny, Nz. 
The other parameters are the same as the Table U 


Lx 

4.373 X 10“® 

[m] 

Ly 

4.373 X 10“® 

[m] 

Ai 

2.000 X 10^ 

[ms-2] 

A 2 

6.000 X 10^ 

[ms-2] 

Nx X Ny X Nz 

64 X 64 X 64 



case of the square pattern. However, the information becomes crucial in the case of more complex patterns 
as we will see later. The resultant grid on the Fourier space is shown in Fig.[6ja). 

We start the simulation with zero velocity everywhere and the perturbed flat interface. The perturbation 
of the interface is given in terms of the Fourier modes = 0) for the wavenumber range 0 < \kx\ < 

0.5max(fca;) = ttNx/{2Lx) and 0 < \ky\ < 0.5max(fcy) = 'KNy/{2Ly). The real and imaginary parts of 
C(fc) in the range are set by independently and identically distributed random variables with a uniform 
distribution between —1/2 and 1/2. The zero Fourier mode is set (/(fc = 0) = 0.2Lz. We lastly transform 
/(fc) in the physical space and multiply the perturbation by an arbitrary factor so that max(|(/(a:) — 0 . 2 ^ 2 !) = 
5.0 X lO-^L^. 

We use here both the original level-set method and the particle level-set method for comparison. 

As shown in Fig. [5J indeed a square pattern is obtained in the simulations. With both the original 
level-set method and the particle level-set method, we start to recognize the square pattern around t ^ Ty. 
In spite of the same appearance of the pattern, the long-time behaviors of the two level-set methods are 
different. With the particle level-set method, the temporal variation of the interface elevation at a point 
(a;, y) = {0.5Lx, 0.78Ly) reaches a steady state around t ^ 38T„ as seen in Fig. [3] (solid line). In contrast, 
with the original level-set method it does not reach a steady state but keeps increasing as depicted with the 
dotted line in Fig. [31 However the square pattern is not destroyed by the unsteadiness up to t = 45r„, at 
which we end the simulation. 

To characterize the difference between the original and particle level-set methods we here introduce two 
time scales: first pattern recognition time and saturation time. The former is the time we first recognize 
the expected pattern, which is T„ for the square pattern case. The latter is the time needed to reach the 
steady state, which is 38T.u for the particle level-set method. Although these time scales are determined 
subjectively and are dependent on the initial condition, they play a useful role in comparison between the 
two level-set methods as we will discuss later. 

Our next target is the hexagonal pattern observed in the experiments^di. The physical parameters of the 
simulation are listed in Table uni As in the square pattern case, we set the size of the horizontal domain to 
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Figure 2. The interface profile for the square-pattern regime at t = 44.5T„. The interface is colored according 
to its height: the red corresponds to higher region and the blue corresponds to lower region. Each horizontal 
length of the domain displayed here is three times that of the calculation domain. The aspect ratio of this figure is 
{3Ly)/{3Lx) = 1.000. Here we show the result with the particle level-set method. 



Figure 3. Temporal variation of the interface height at position {x, y) = {0.5Lx, 0.78Ly) for the square pattern 
calculated with the particle level-set method (solid line) and with the original level-set method (dotted line). 


the minimal size containing one hexagon. Specifically, with the resonant wavevector k[ shown in Fig. |6Kb), 
the lengths are = 2TT/k[y. and Ly = 2Trlk[y. The aspect ratio is Ly/L^ = 0.5774 r; 1/-\/3- The numbers 
of grid points used per wavelength are {Xi/Lx)Nx = 0.5Nx and {Xi/Ly)Ny = 0.8660A^y ~ \/i/2Ny. We use 
the same initial condition as the square pattern case. 

The hexagonal pattern is reproduced with both level-set methods. The result with the particle level-set 
method is shown in Fig. 21 Despite the pattern being the same, the first pattern recognition time is different 
between the two level-set method: QTy for the original level-set method and IQTu for the particle level-set 
method. The saturation time is 24T„ with the particle level-set method as shown in Fig. [5] In contrast, 
saturation does not occur with the original level-set method during our simulations of length t = 40T„. The 
hexagonal shape of the pattern is maintained in spite of the unsteadiness. 
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Table III. Parameter values for the hexagonal pattern. The other parameters are the same as the Table H] 


Lx 

Ly 

Ai 

A2 

X Ny X 

1.184 X 10“^ 

6.837 X 10“® 

3.200 X 10^ 

3.000 X 10^ 

112 X 64 X 64 

[m] 

[m] 

[ms-^] 

[ms-2] 

Table IV. Parameter values for the 2k rhomboidal states. These parameters (except for Lx,Ly, 
Nz) are identical with the experiment by Arbell et al^ 

Lz,0, Nj:, Ny and 

Lx 

1.446 X 10"2 

[m] 

Ly 

5.234 X 10"® 

[m] 

Lz 

1.000 X 10"^ 

[m] 

Pt 

1.293 

[kgm-®] 

Pb 

9.500 X 10^ 

[kgm-®] 

Vt 

1.822 X 10"® 

[kgm-®s-®] 

Pb 

2.185 X 10"^ 

[kgm-®s-®] 

Ai 

2.372 X 10® 

[ms-2] 

A2 

4.925 X 10® 

[ms 2] 

UJl = 2 (jJq 

3.141 X 10^ 

[s-®] 

0 J 2 ~ 3(Uo 

4.712 X 10^ 

[s-^] 

e 

0 

[rad] 

a 

2.150 X 10"^ 

[kgm-®] 

g 

Nz: X Ny X Nz 

9.807 

96 X 32 X 64 

[ms-2] 

bottom-fluid depth 

2.00 X 10"® 

[m] 


The above results on the square and hexagonal patterns suggest that the original level-set method is not 
a suitable interface-tracking scheme for Faraday waves. Although the patterns initially emerged with the 
original level-set method are consistent with the experiment, it is seen that the temporal variation of the 
interface height does not reach a steady state. This unsteadiness in the long run may change the correctly 
selected pattern initially into a different shape with the original level-set method. In the simulation of the 
rhomboidal pattern, the deficiency of the original level-set method appears more seriously as we see in the 
next section. 


C. Rhomboidal states 

The next pattern we seek to simulate is called the 2k rhomboidal state observed in the experiment by 
Arbell et al^. The pattern is observed around the bicritical point which appears as the sharp tip in Fig. [TJ 
There are two linearly unstable wavenumbers ki and ^ 2 , hence the name 2k rhomboidal states. As a result 
of the nonlinear interaction among the resonant modes a simple resonance relation appears: + ^2 = ki, 

as shown in Fig. Hie). This rhomboidal pattern involves two circles in the wavenumber space, which is a 
notable difference from the square and hexagonal patterns. 

The experiments on the rhomboid patterns were reported in the two references^id^. There is a slight 
difference in the experimental settings between the references. We succeed in simulating the rhomboidal 
patterns with the same parameters for each of the two references. However here we present only the 
result corresponding to one of the references^ since it contains a detailed analysis of the pattern along with 
a photograph of the rhomboidal pattern. Note that for the square and hexagonal patterns we use the 
parameters of referencei^. The numerical parameters are listed in Table HVl We set the size of the horizontal 
domain again to be minimized containing one rhomboid, namely = 27r/(fci/2), = 27r/(A:2 sinyj) where 

(fi = 70.16 is the angle between the vectors fci and k 2 shown in Fig. [UJc). It is calculated from the relation 
2k2 cos(p = ki- The aspect ratio is thus Ly/L^ = 0.3620. The numbers of grid points used per wavelength are 
{2TT/ki){N,/L,) = 0.5N,, {2TT/ki){Ny/Ly) = 1.387Ny, {2TT/k2){N,/L,) = 0M08N, and {2TT/k2){Ny/Ly) = 
0.9415fVy. The initial condition is set in the same way as the cases of the square and hexagonal patterns. 
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Figure 4. Interface profile for the hexagonal pattern at t = 27.65r„(above), 28.22T.„(bellow). The interface is colored 
according to its height: the red corresponds to higher region and the blue corresponds to lower region. The dimension 
of the domain displayed here is SZ/a, x 3Ly x Lz- The aspect ratio is LyjLx — 0.5774 « l/v^. Here we show the 
result with the particle level-set method. 


With the original level-set method, we do not obtain the rhomboidal state. On the other hand, with 
the particle level-set method, we obtain the state as a steady state as shown in Fig. |8l The first pattern 
recognition time of the rhomboid with the particle level-set method is t ^ 29Ty and the saturation time is 
the same t ^ 29T„ as depicted in Fig. 0 The first pattern recognition time is much longer than those of the 
square and hexagonal patterns. We consider that the nonlinear interaction among the resonant modes on 
the two circles in Fig. |6l(c) takes a longer time in order to reach a constant oscillation amplitude. 

Figure [5] shows the temporal evolution of the Fourier amplitudes of the interface height for the three 
resonant modes. The circle symbols represent the simulation results and the solid line is the evolution 
calculated with the Floquet coefficients obtained in the linear stability analysis^^. The evolution of the 
nonlinear rhomboidal modes (circle symbols in Fig. |9]) is quite close to that of the linear results, which 
indicates that the nonlinear effect in the temporal evolution of the pattern is weak. 
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Figure 5. Temporal variation of the interface height at position at (a;, y) = (O.SLx, 0.5Ly) for the hexagonal pattern 
calculated with the particle level-set method (solid line) and with the original level-set method (dotted line). 


Now we compare the simulation results with the weakly nonlinear analysis of the rhomboidal states by 
Porter et Their analysis for the first time explains with an elegant broken-symmetry argument why 

the rhomboidal pattern appears. In deriving their amplitude equations up to third order in the amplitude, 
they assume that the rhomboidal state is close to the bicritical point and that the damping parameter 7 
is small. Accordingly they expand the coefficients in the amplitude equations in powers of the vibration 
amplitudes A* = (Ai — Aic)/Aic, A 2 = {A 2 — A 2 c)fA 2 c and the damping parameter y/wo- The resulting 
coefficients of the quadratic term of the amplitude, their signs and dependence on 7 , explain the rhomboidal 
pattern selection for certain frequency ratios m : n. However, as they discussed, it is not clear that the 
damping parameter y/wo is small enough in the experiments^. 

To test the assumption, we measure the damping parameter from our simulation data. Before doing this 
we calculate it with dimensional analysis: the damping parameter of the bottom fluid can be estimated as 
'y{k)/u!o = 2r]b/{pbk‘^uJo) with the critical wavenumber k. This gives ■j{ki)/uJo = 0.219 and ‘j{k 2 )/uJo = 0.476, 
where the critical wavenumbers ki and ^2 are determined for the critical vibration amplitudes Aic = 21.9 and 
A 2 C = 44.8. These dimensional values can differ in orders of magnitudes from the actual damping parameter. 
In our nonlinear simulation of the rhomboidal pattern, we set the normalized vibration amplitudes A* = 
(23.7 — 21.9)/21.9 = 0.0831 and A 2 = (49.3 — 44.8)/44.8 = 0.100, which justifies the expansion in terms of 
Al and A 2 in the coefficients of the weakly nonlinear analysis. In order to estimate the damping parameter 
in our nonlinear simulation, we used a method similar to that used in the experiment^: we take the 
snapshot at t = 141.4Ti, from the rhomboidal pattern simulation as an initial condition; we then start 
the simulation without the vibration forcing and measure how the interface elevation decays in time. The 
temporal interface behaviors on the line at xj = 0 are shown in Fig. HOI The envelope in the figure gives 
')Ty ~ 1.23. Consequently, the damping parameter y/wo is 0.196. Although this is smaller than unity, it 
may not be small enough to ignore its higher order terms. 

We next try to obtain the slowly varying amplitudes from the fully nonlinear evolution of the resonant 
modes shown in Fig.|9l For example, ImC(fe 2 ) divided by the sub-harmonic oscillation C sin[27r/(w2/2) t-l-0], 
where 0 J 2 = 3uJo, C is a suitable factor and 0 is a suitable phase, should give a slowly evolving function. 
We divided the resonant mode (symbols in Fig. [5]) by the sub-harmonic oscillation. However the calculated 
function is not slowly varying in time. Moreover, we divided the nonlinear data (symbols) by the linear 
Floquet-mode data (lines) in Fig. |9l The calculated function is not slowly varying either. Nevertheless 
we look at the phase-space orbit formed by the three variables in Fig. We do not find a characteristic 
structure often associated with the solutions of the normal-form equations corresponding to the rhomboidal 
structure. Hence we are not able to compare our data with the weakly nonlinear analysis in this respect. 
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Figure 6. Nonlinear resonant wavevectors (a) square pattern: |fei| = |fci| = 1.436 x 10® (b) hexagonal 

pattern: |fci| = = |fci| = 1.061 x 10® the angle between k'^ and the axis is 60°. (c) rhomboidal 

pattern:|fei| = 8.653 x 10® |fc 2 | = 1^21 = 1.275 x 10® the angle ^ between fc 2 and the axis is 70.16°. 

The grid represents the minimal discretization of the Fourier space for each pattern. 



Figure 7. Temporal variation of the interface height at position at (x, y) — {0.3Lx, 0.25Ly) for the rhomboidal state 
calculated with the particle level-set method (solid line) and with the original level-set method (dotted line). 


D. Comparison between original and particle level-set method 


We observe that the original level-set method and the particle level-set method yield qualitatively dif¬ 
ferent results. With the original level-set method, the square and hexagonal patterns are observed but do 
not become constant-amplitude oscillations. The rhomboidal state, which is here the main target, is not 













































































Numerical simulation of Faraday waves oscillated by two-frequency forcing 


16 



Figure 8. Interface profile for the rhomboidal state at t = 40.38T„. The interface is colored according to its height: 
the red corresponds to higher region and the blue corresponds to lower region. The dimension of the domain displayed 
here is 2Lx x ALy x L^. The aspect ratio of the displayed domain is ‘iLy/{2Lx) = 0.7239. Here we show the result 
with the particle level-set method. 




Figure 9. Temporal evolution of the resonance amplitudes of the interface height C, for the wavevectors fe 2 , fci 
(see Fig. [HKc)) in the rhomboidal state. Here we show only the dominant parts for each wavevector. The symbols 
are data of the simulation with the particle level-set method. The solid lines are time evolution of the neutral stable 
modes calculated with the ten Floquet coefficients in the linear stability analysis. 


observed. On the other hand, in our simulation with the particle level-set method, all three patterns are 
observed and become constant-amplitude oscillations in agreement with the experiments. This difference is 
due to the well-known problem of the original level-set method, which we discuss here. 

In order to clarify the difference between the two level-set methods, we look at how well the volume of 
the lower fluid is conserved during the time evolution. The volume of lower fluid is calculated with H((f>), 
Eq. (jlOp . as V{t) = f H{(j){x,t))dx. The variations of the volume for the hexagonal and rhomboidal cases 
are shown in Figs. [TTJ [T^ with the numerical parameters listed in Tables MM 

As shown in Figs. [TT] and m the volume increases with the original level-set method, instead of being 
conserved. This non-conserving property of the original level-set method is well known^^— . This explains 
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Figure 10. Measurement of the damping parameter 7 . The curves are temporal variations of the interface height 
at various points on the line 1 = 0 without the vibration forcing. The damping parameter 7 /tJo can be estimated 
from the envelope 0.28exp(—1.23t/T„). The interface height at the quiescent state (where the Faraday waves decay 
completely) is denoted as QeqjLz — 0.197, which is slightly different from the initial interface height 0.201/z without 
the random perturbation. 




1.12 

1.10 

1.08 

1.06 

1.04 

1.02 

1.00 





— particle level set method 

- - original level set method 


0 5 10 15 20 25 30 35 40 

i/7; 


Figure 11. Comparison of variation of the bottom fluid volume between the original level-set method and the particle 
level-set method in the case of the hexagonal pattern. 


why the interface height does not reach constant-amplitude oscillations with the original level-set method 
for the square and hexagonal patterns. Concerning the rhomboidal pattern, the original level-set method 
fails to exhibit the pattern. But with the particle level-set we start to recognize rhomboids at t = 29r„ (first 
recognition time). At this time it is seen from Fig. [T^lthat the volume in the simulation with the original 
level-set method increases by 10%. In other words, a long time is needed for the nonlinear interaction to 
form the resonant modes for the rhomboidal pattern. During this time, the error of the simulation with the 
original level-set method, the increase of the bottom-fluid volume, becomes so significant that the rhomboidal 
pattern is not observed. Therefore we conclude that the particle level-set method is more suitable than the 
original level-set method to reproduce complex patterns such as the rhomboidal pattern, which require a 
long time for selection 
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Figure 12. Same as Fig. 1111 but for the case of the rhomboidal states. 


IV. SUMMARY AND DISCUSSION 

Motivated by the recent experiments of Faraday waves with two or more frequency forcings exhibiting 
even richer patterns than the single frequency case, we have conducted a numerical simulation of the two- 
frequency Faraday waves, specifically targeting the rhomboidal pattern. 

We first validated our numerical simulation with the linear stability analysis of the two-frequency Faraday 
wavesi^. The two simple patterns, the square and hexagonal patterns, in the nonlinear regime were simulated 
with the same physical parameters as the experiments^. In particular, the simulation using the particle 
level-set method in the minimal computational domain reproduced the two patterns in agreement with the 
experiment. Employing the particle level-set method, we finally reproduced numerically the 2k rhomboidal 
states, the most complex pattern in this numerical study, with fluid properties identical to those of the 
experiments. We next checked whether the rhomboid obtained in our simulation satisfies the assumption 
made in the weakly nonlinear analysis for the rhomboidal patteroSSiS^. Specifically, the assumption concerns 
the smallness of the damping parameter and the vibration amplitudes. We found that the damping parameter 
of the rhomboidal pattern in our simulation is marginally small. Further comparison with the weakly 
nonlinear analysis is difficult. 

In these simulations, we used two level-set methods: the original level-set method and the particle level- 
set method. The interface motion of the Faraday waves appears quite modest in the sense that it is not 
usually considered as a typical target of the interface-tracking schemes. One may think that any modern 
scheme is capable of simulating Faraday waves. However, due to the well-known problem of the original 
level-set method^, we failed to simulate the square and hexagonal patterns as steady states and to reproduce 
the rhomboidal pattern at all. Thus the Faraday wave problem requires an accurate scheme tracking of the 
interface such as the particle level-set method. One reason for this is that we need to simulate the system for 
a long time if we start with a random initial condition. We believe that, in developing a new implementation 
of the interface-tracking scheme, the Faraday wave problem can be a benchmark problem in addition to a 
physical phenomenon. In the linear regime quantitative comparison can be made as demonstrated in the 
simulation by Perinet et alr^. In the nonlinear regime qualitative comparison can be made (whether or not 
the right pattern emerges if we choose parameters for a certain pattern observed in experiments). 

We carried out simulations on the minimal calculation domains to reproduce the three patterns. Its effect 
was studied for the rhomboidal case in the following way. Simulations were run in domains which were twice 
as large with twice as many points, thus keeping the density of numerical grid points constant. Accordingly 
we have the same grid spacings in the physical space, L^/N^ and Ly/Ny. (for the z direction we keep the 
same values for Lz and Nz). The rhomboidal pattern is observed with this setting with the same first pattern 
recognition time and the saturation time. Hence it is unlikely that the minimal domain setting affects the 
pattern selection numerically. We also checked whether the number of grid points in the vertical direction 
Nz is sufficient or not by doubling Nz but retaining the other parameters as in Table IIVI The result does 
not change. 

As we mentioned briefly in the Introduction, many other patterns are observed in the experiments of the 








Numerical simulation of Faraday waves oscillated by two-frequency forcing 


19 


two-frequency forced Faraday waves. In fact, our initial goal was to reproduce not only the rhomboidal 
pattern but also the hexagonal based oscillon (HBO) (also known as double hexagonal superlattice (DHS)), 
the spatially subharmonic superlattice pattern (SSS) and the oscillon observed experimentally by Arbell et 
So far, we have not been able to reproduce those patterns perhaps due to our strategy to use the 
minimal calculation domain including one pattern. Setting the minimal domain corresponds in terms of the 
Fourier space to maximizing the grid spacings in the kx and ky directions to include the resonant modes 
with discretized points. For these patterns we failed to simulate; in fact, it is not clear how to set a minimal 
domain even with the knowledge of the selected resonant modes available from the experiments. 

Now we take the HBO pattern as an example and discuss the difficulty of setting the minimal domain. 
In the linear analysis of the HBO case, two different wavenumbers simultaneously become unstable. Hence, 
as in the case of the rhomboidal pattern shown in Fig. [51 two circles can be important. However according 
to the experiment the resonant modes lie on only one of the two, which we call the resonant circle; we call 
the other circle the non-resonant circle. As a first trial we took the minimal calculation domain to resolve 
only these resonant modes on the resonant circle without including modes on the non-resonant circle. With 
this minimal domain and the particle level-set method, we did not obtain the HBO pattern at all starting 
from the same initial condition as in Sec. IHIl We speculate that, in the course of establishing the resonant 
modes, the modes on the non-resonant circle are important in the pattern selection and hence should be 
taken into account properly in the simulation. Of course, if we could enlarge the calculation domain and 
increase the number of grid points in the physical space in order to take a large number of mesh points 
near the non-resonant circle in the Fourier space, this problem might be overcome. Even though we have 
doubled Lx, Ly, Nx and Ny, in order to take smaller grid spacing in the Fourier space, the HBO pattern 
did not emerge. A finer grid would make the cost of computation prohibitively high (notice that a long time 
simulation is also needed here). 

We also ran the simulation of the SSS but failed possibly for the same reason. The resonant modes of 
the SSS observed experimentally lie either on a circle whose wavenumber (radius) is linearly stable or one 
of the two circles determined by the linear stability analysis. For both the HBO and the SSS cases, we 
checked that the volume is conserved to the same degree as it is in the rhomboidal case with the particle 
level-set method. Regarding the oscillon, the structure of the resonant modes in the Fourier space is not 
clarified experimentally, implying that we do not have any guidance on the discretization of the Fourier 
space. Perhaps, guessing from the physical-space appearance of the oscillon, the number of excited Fourier 
modes is very large compared with other patterns. To circumvent this sort of difficulty, a completely different 
numerical scheme with Chebychev polynomials for capturing a localized structure is proposed by Lloyd et 
ali^S, which may be worth exploring. Moreover the oscillon’s metastabilityi^ may make simulation even 
more challenging. Our future work is an approach relying on computing power in which we take as high a 
resolution as possible to reproduce complex patterns like the HBO, the SSS and the oscillon. We believe that, 
if such a simulation succeeds, it would provide knowledge about the role of the modes on the non-resonant 
circles in pattern selection. 
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